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FaCE is a self contained programme, with namelist input, that solves the three body Faddeev 

equations. It enables the inclusion of excitation of one of the three bodies, whilst the other two 

remain inert. It is particularly useful for obtaining the binding energies and bound state structure 

compositions of light exotic nuclei treated as three-body systems, given the three effective two body 

• • ' ■ interactions. A large variety of forms for these interactions may be defined, and supersymmetric 

^i ' transformations of these potentials may be calculated whenever two body states need to be removed 

^^ I due to Pauli blocking. 

^ ■ PACS numbers: 11.80, 21.10, 21.45+v, 21.60.Gx 

C ^ Keywords: three body problem, core excitation, exotic nuclei, bound states, Faddeev equations, hyperspher- 

3 , ical harmonics 

CN : PROGRAM SUMMARY 

, I Title of program: FaCE (Faddeev with Core Excitation) 
^ Computers: The code is designed to run on any unix/linux workstation or PC. 

Operating systems: Linux or UNIX 



f^ Program language used: Fortran-90 

Numerical libraries used: Source code for 6 routines from the NAG and BLAS libraries is included to enable indepen- 
dent compilation. 



o 

f^ Memory required to execute with typical data: 9 Mbytes of RAM memory and 12 MB of hard disk space. 

r^ No. of bits in a word: 32 or 64 

I , No. of lines in distributed program, including test data, outputs, etc.: 13944 

O ' Distribution format: compressed tar file 

Keywords: three body problem, core excitation, exotic nuclei, bound states. 

Nature of physical problem: The program calculates eigenenergies and eigenstates for the three body problem by 
solving the Faddeev equations. 

^\ Method of solution: Given the two body effective potentials it performs the supersymmetric transformation in case 

where there are forbidden states to be removed. The three body wavefunction is expanded in hyperspherical coor- 
dinates, the hyper-angular part is a series of jacobi polynomials and the hyper-radial part is written in terms of a 
laguerre basis. Within this basis the three body matrix elements are calculated and the full three body Hamiltonian 
matrix is completed. The diagonalization process is performed after various reductions (isospin, orthonormal and 
Feshbach) to determine the energies. Finally the three body wavefunction is reconstructed and other bound state 
observables are calculated. 

Typical running time: 6 sec on a 1.7 GHz Intel P4-processor machine. 
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LONG WRITE-UP 



I. INTRODUCTION 



Radioactive nuclear beams have allowed the exploration of the nuclear driplines (proton rich and neutron rich), 
and unveiled exotic phenomena. Many of the properties of light exotic nuclei have been well described within few 
body formalisms: one neutron halos such as ^^Be and "'^^C have been described with two body (core+N) models 'l','^; 
Borromean systems such as ^He and ^^Li have been well modelled as three body (core+N+N) systems j|, ,4j; and ^He 
has been successfully accounted for within a five body picture 5j . Although in the early days these models assumed 
all participants were inert, and concentrated on treating the few body dynamics exactly, the advantage of retaining 
degrees of freedom of the core was soon realised |g|. Few body wavefunctions, solutions to the few body Hamiltonian 
with core excitation, should contain the main components of any microscopic calculation, with the advantage of its 
simplicity. 

Few body models have become extremely useful in the field of light radioactive nuclei, not only from the structure 
perspective, but mainly for the purpose of reaction modelling Q. Some important consequences were found when 
extracting radii from reaction cross sections Q , when analysing transfer reactions for extracting spectroscopic factors 
I^J or when studying elastic and inelastic scattering [l^l . Many of the features contained in the few body structure 
models are essential for a good description of the reaction process. 

In this paper, we present a self contained program that provides a solution to a general three body problem where 
one of the clusters is allowed to excite. In Section II a brief overview of the construction of the three body basis is 
presented. In section III the matrix elements required for the standard interaction are given. Section IV discusses 
matrix reduction methods (useful for big calculations) and then the diagonalisation procedure. Section VI contains a 
summary of the observables that are calculated in FaCE. In Section VII specific comments on the program and the 
input manual is provided. Finally in Section VIII we illustrate the use of FaCE with three physical examples. 



II. THE THREE BODY BASIS 

Our intention was to develop a general tool to handle the bound state properties of a nucleus well described as 
a three body system i+j+k where one of the particles is allowed to excite. FaCE is based on solving the Faddeev 
equations [ll| with a hyperspherical formulation of the general three body problem [3, U2, UM ■ 

The Faddeev equations define three component wave functions "^(^^ ^ such that the full three-body wavefunction 
is xji-'^^ — ^{^{xi^yi) + ^2*^(2^2,2/2) + ^3*^(a^3>y3)- Here, the components ^i are functions of their own 'natural' 
Jacobi coordinate pairs i (as in Fig. ^, and are solutions of the Faddeev coupled equations: 



{Ti + h + Vi- 


-E)^'i^' = 


-T^i(*f' + *fO 


{T2 + h + V2- 


-E)^i^^ = 


-T^2(1'3*^ + *5'*0 


{Ts + h + V3- 


-~E)^i^^ = 


-Vsi^i^^ + *f 



(1) 

These equations contain h — ^^hi, the sum of the intrinsic Hamiltonians of each particle hi, the relative kinetic 
energies in each coordinate set T, = T^i + Tyi and the two body interactions between the corresponding pair Vi = 
^jki^jk) (both the Coulomb and nuclear interactions). The indexes i,j,k run through (1,2,3) in circular order. 

The distances between each pair of particles fjk, and the distance between the centre of mass of the pair and the 
corresponding third particle (represented in Fig. (^ by the thin lines), can be expressed in terms of the Jacobian 
coordinates {xi,yi) where Xi = ^3^ f^-fe and -ffi = yZ^^fQi,),,: 

^j/c = ^i ~ He and fQ-j,),; = t^j - {Ajfj + Akrk)/{Aj + Ak) ■ (2) 

Note that the reduced masses are defined by Ajk — ^',X ^^'^ ^O'fe)* ~ a+a+a' ^^^^ hj^k G (1,2,3) where 
^^ = Hii with 771 = 1 a.m.u. and m^ the mass of particle i in a.m.u. In FaCE we will use X to refer to the pair 
{xi,yi), Y to refer to (2:2,2/2) and T to refer to (x^^y^). 

FaCE allows the user to include core excitation of one of the particles. Let us assume one particle c S (1, 2, 3) has 
low lying excited states strongly coupled to the ground state, and which are likely to have important roles in the three 
body system. Particle c is then treated as the core, and its internal coordinates ^c must be added to the set of Jacobi 









FIG. 1: Three sets of Jacobi Coordinates used in the Faddeev Formahsm. 



coordinates to define the full quantum state of the system. The intrinsic Hamiltonian of the core determines a set of 
eigenstates (f>s^ and eigenvalues e^^ , 



The model then expands the total wavefunction of the system in terms of these 
degrees of freedom from the Jacobi coordinates in each Faddeev component: 



(3) 
states, and factorizes the core 
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{xi,yi,ic)^^ ipsAic) ^sAxi,yi) 



with i = 1,2,3. Here V'sc contains the radial, angular and spin of the remaining two particles relative to the chosen 
core. This model is advantageous if only a small number of core states ^s^ is required to describe the system accurately, 
which is normally true for systems close to the driplines. 

FaCE uses the hyperspherical method to convert two-dimensional partial differential equations into a set of coupled 
one-dimensional equations. The Jacobi coordinates (xi, yi) are transformed into the hyperspherical coordinates (hyper- 
radius Pi and hyper- angle Oi) defined as 



= E^» 



and 9i — 



arctan( — ) 
yi 



(4) 



The hyper-radius is invariant under translations, rotations and {i,j) permutations, and is directly related to the overall 
size of the nucleus whereas the hyper-angle contains radial correlations and is related to the relative magnitude of 
the two Jacobi coordinates. The hyper-radius is the same for all i=l,2,3, this being a basic advantage offered by the 
hyperspherical coordinate system, while the hyper-angle 9i = arctan( — ) is different for the various X, Y, T bases. 

The transformation from a Jacobian coordinate set to a hyperspherical coordinate system does not affect the 
angular and spin variables, nor the degrees of freedom of the core. Isospin dependence is not explicitly introduced 
since typically the interactions to be used have a fixed isospin. 

For a given Jacobi set [xi,yi), we need to define and couple together the associated orbital angular momenta 
{Ixi, lyi), as well as the states and the spins of the three particles Si, Sj, Sk, as shown in Fig. |21 The core particle will 
have in FaCE an index for its different excited states, but for the presentation below we assume that the spin value 






FIG. 2: Jacobi coordinates and the notation for the corresponding angular momenta. 

Sc suffices to identify its state. A partial wave decomposition for each Faddeev component uses the following the 
coupling order 
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JAI 



J2^a■ix^,y^) \i : a^Y 



(5) 



with the abbreviation Ui = {{Ixi, lyi)Li, (sj, Sk)Sxi}Ji] Si for the quantum numbers of each component i. One of the 
Si, Sj, Sk in ai will identify the excitable core state Sc- 

The two-dimension radial wavefunction ip^^ {xi, yi) is next expanded in the hyperspherical variables. The separation 
between hyper-angle and hyper-radial dependence of the wavefunction makes use of the fact that the hyper-angle 
functions, eigensolutions of the hyper-angular (centrifugal) part of the three-body kinetic energy operator [I3, are 
explicitly defined in terms of the Jacobi polynomials: 



C':(a;^,y» 



^Y^^kM vt'^i) 



Ki 



with LpY^^Oi) = iY^*'"' (sine*)'" (cos 6')''" Pl^;+^l^'^y^+^l^(cos26,) 



(6) 



(7) 



Here P„i is the Jacobi polynomial, NKi is a normalisation coefficient and Ki is the hyper- angular-momentum directly 
related to the order of the corresponding Jacobi polynomial Ki — Ixi + lyi + 2ni (ni=0,l,2,...). In order to sim- 
plify the notation we will omit whenever possible the total angular momentum and projection labels JM from the 
wavefunctions. 

Introducing this expansion in the Faddeev Equations, and performing the hyper-angular integration, one obtains a 
set of coupled equations for the wave functions x^a k (p) °^ ^Q- ® 



(Tp + Lk^ ip) - i?)xL.K. (P) + E Kk,,c.,k, (p)xi,K, (P) = 0, 

jaj Kj 



(8) 



^ ''■ and the centrifugal potential is LK^ip) = h'^{Ki + 3/2){Ki + b/2)/{2mp'^). The couplings are 



where Tr, — „ . ■ 

f^ 2m dp 

the hyper-angular integrations of the two-body interaction V^''^. a k (p) — < 'fix- ^^ i^j)\^ij\VK- "'(^j) >• In FaCE, 
these hyper-angular integrations are performed using Gauss- Jacobi quadrature on a grid with Njac points (defined in 
namelist grids in the manual). Gauss- Jacobi quadrature points are evenly spaced in hyper-angle. 

In order to solve these coupled equations, the hyper-radial behaviour is expanded in terms of orthonormal basis 
functions 



i?„(p) = p'/^Po^[nl/{n + 5)f/'Ll{z) exp(-z/2) , 
where z = p/po with scaling radius po £ind L'^{z) is an associated Laguerre polynomial, as 



(9) 



(10) 



The potential matrix element integrals of the i?„(p) functions are calculated using Gauss-Laguerre quadrature with 
Niag points, which must be greater than the number Nf, of basis polynomials. Nf, is set from the namelist solve, 
and Niag is set from the namelist grids, while the quadrature points and weights are determined through finding 
numerically the roots of i^ (z) = 0. 

The kinetic energy matrix elements, including the centrifugal barrier, are 



{R,^{p)\TK,\Rn'{p)) = TT 

Im 



n< , KdK^ + A) 
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{5(n> - n< + 1) + n> + n< + 1} 



where n< = min(n,n') and n> = max(n, n'). 

After introducing the hyperspherical expansions Eqs. (|6ll0f) into the Faddeev coupled equations Eq. Q, one arrives 
at a set of simultaneous linear equations 



Hsi^Esi 



(11) 



for the coefhcients a = {a}. We shall only be calculating bound or pseudo-bound states, for which the wave functions 
x{p) vanish at both p = and p — > oo. This is guaranteed by those same properties of the basis functions of Eq. ^. 



III. THE THREE BODY MATRIX ELEMENTS 



The complete wave function solution for a given J (which will henceforth often be omitted) is 



* = I]I]C(2^»'yO I*: "^) 



(12) 



i—\ cti 



where there is an implicit sum over hyper-moment K due to the expansion of V'^. (x^, y^) as in Eq. ©. The Hamiltonian 
matrix will therefore require overlap integrals of the potentials between pairs of the overcomplete basis set {|i : ai)}. 
We will need transformation matrices for the rotations \k : Uk) -^ \i '■ OLi) clockwise, and \i : a^) ^- \j : aj) 
anticlockwise, between the three Faddeev components. Considering (i,j,k) the circular order, the expressions that 
allow the transformation (which conserves total angular moment L and hypermomcnt K) between Faddeev components 
in both directions are 



^2{J-S^k-Sk)+St+S^ 
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(13) 



(14) 



where RR{lxi,lyi',lxjTlyj\L) are the Raynal-Revai coefficients [lj|. The two above equations introduce two kinds of 
norm matrices N^^ and iV*'"' such that \i) — N^-'\j) and \i) = N'^^\k). That is, the matrix elements of N^^ are the 
basis-state overlaps 



iV«J 



{j : aj\i : ai) 



so that mN^'' = N'''. 

There are several kinds of matrix elements needed for the matrix H in the Faddeev equations of Eq. Ijllll . The 
general potential we are considering is: 



Vjk = V,ix,) = V^ix,) + SO Vr{x,) + Q V:^{x,) + T V/ {x,) + SS V^^x,) 



(15) 



where V^ stands for the central interaction, SO and V^° are the spin-orbit operator and the spin-orbit radial form- 
factor respectively, Q and V^ (xi) are the tensor operator and radial shape for the multipoles of the deformed potential, 
T and Vi^{xi) stand for the standard tensor operator |15|| and radial dependence for the tensor NN interaction, and 



finally SS and Vf^ are the spin-spin operator and the corresponding form-factor. All these are included in FaCE. 
The parameters for the corresponding radial form factors are defined in namelist poten (see manual for details). 

The matrix elements oiVi{xi) are preferentially calculated between basis states of the same i component. The norm 
matrix elements N^^ allow us to express general potential matrix elements in mixed representations, in terms of the 
preferential Faddeev representation. For example: 

{j : a'^m\z : a,) - E ^1^ (* ■ "^H^ I* ■ "^) ' (16) 

a'. 

or, for cases when the potential is most easily presented in one particular Jacobi coordinate set: 

(/ : a'^mWj : a,) = E N% {i : a'M¥ ■ c^^N^., ■ (17) 

— I j 

In this section we first consider the angular and spin matrix elements, and these will be later multiplied by numerical 
integrals over the hyper- angle 6i to obtain V^j^, ^ k (p) ^^ ^^ ^'i- ® ^^'^ ^^^ hyper-radius p after the expansion in 
Eq. (|10l) . We will use Wu' to denote the matrix elements over the angular momentum basis states. We use: 

{i : a'M¥ ■■ «.) = Wt,, V^ix.) + W^ Vr{x^) + W^^, V.'^ix,) + Wj, V^ [x,) + WfJ^V^^^x,) . (18) 

The potential matrix clement for the central part is diagonal in all angular and spin variables. Next, let us 
consider the spin-orbit part. As all three particles may have spin, we have introduced the general spin operator 
Si = rySj+FjfcSfe, where Fy and F^fc select which of the spins are to be dynamically coupled, and with which relative 
strength. The matrix elements for this operator are 

L^ h 1 1 f L\ U 1 

"ixi t^xi "^i ] \ '2:1 I'xi ''yi 

< Qi <? . 1 1 1 

.(19) 



>^xi >^^-; ^i ^xi ''-^■i ^i. 



ry(-i)''"s,Ys,(s, + i)|^^' "^r ,^^} + r.fc(-i)^^-4V5fe(sfc + i)| 



'-'7;?: ^x'i' ^ 
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Sk S'fc s'- 



Typically the interaction between a deformed excitable nucleus and another particle is expanded in multipoles. 
The essential angular momentum operator is a tensor interaction of the type Cq{Ix) ■ Cq{s). Depending on the 
Faddeev component and the subscript that specifies the deformed nucleus, all forms of the operator are needed: a) 
Cqllxi) ■ CQ{si), b) CQ{lxi) ■ Cq{sj) and c) Cqllxi) ■ CQ{sk)- We next present the results for these matrix elements: 

a) {i : a[\\CQ{li) ■ Cq(s,)||z : a,) = <5,,,/,, ,,5jj,5i,^,^^55^5„(-l)^'+''''+'-+'''+'^'+^''+''"i;j^i^i»4CzS^ 

x{5 5 ?}{ii ill, }{££;?}( I ?'s')<'"«^<*'^ '^»' 

h) {i : a[\\CQ{ix^) ■ Cq{s,)\\i : a,) = 5,;,,(5,. ,,(5jj,(5j,j.5r _i^^(-l)^'+'U+si>+S-+4+4L^i^5^^5„44s/. 

X { £ 1: 1 } { £ £ ^ } ll' ^' 2 } ( '? 2 »• ) (».ii^(%)ii».> : pi) 

c) (z : a[\\CQ{ix,) ■ Cq{su)\\i ■ a^) = <5,-,^5,-,/jj.,5j^j-^r^^,„.(-l)^'+'-+25"+^-+^^i^i,5V5;,44sl 

^ ^ ^^ - ^ ^ ^ t: £ 2 1 1 !r !r ? K 's' ToO Mc,i^.)w) . (22) 



Jxi ^xi "^i J I ''^* 'xi '•yi J I °k *fc ^j 



The matrix elements between different core states {sii \\CQ{si)\\si) depends on the model used. Under the assumption 
of a pure rotational model, these matrix elements are given by: 



{s,\\CQ{k)\\s.) = {-l^-^'s^ [-b Q k) ^ (23) 

where K is the quantum number for the rotational band. If the interaction is £— dependent there is an ambiguity 
on the choice of the radial form factor (which is defined at the multipole expansion of the deformed potential) . The 
parameter Ipot (see manual) controls this choice. 



When considering the spin-spin interactions again there are three possibihties depending on the Faddeev compo- 
nents: a) (s^||si • Sj\\si); b) (s'j||si • Sfc||si) and c) (s^||sj • Sfc||si). After some algebra one can arrive at: 

Jj Ji 1 1 J ^xi ^xi 1 I J Sj S,; 1 



X y^ii^TTiyy^^iirTT)^ £ s^. l^ ;r 4 i M ^ K' ^ = ^''^ 






(26) 

A reahstic NN force contains a tensor interaction of the type T2{sjSk) ■ C2{lxi) which also needs to be considered. 
Below is the expression for these matrix elements after working out the algebra: 



{i : a',\\T2{s,Sk) ■ C2{lx^\\^ ■ a,) = (5,.,^5,.,/,- ..^jj.Jj.j,^;. j^,(-1)3'^'+'--^"25V5;,£^L,C«44-4\/«j(«j + 1) 






X v^HiTTTTI -^; - ;, 1 1 - £ ,f_ I ^ -y ^ -• j I j^ ., 1 } ,,27, 

Few-body models often include effective three-body potentials to describe the influence of dynamics not explicitly 
described by two-body potentials. We have parameterised the simplest diagonal form of such a potential: 

t(i' ■.a',\\%\\i:a,)^5,,,5^,^^Mp) ■ (28) 

IV. PAULI BLOCKING 

Often within a three-body calculation, it is necessary to eliminate the Pauli forbidden two-body bound states before 
diagonalisation. This may be accomplished by several methods [la. ITi^j: by projection operators inserted in three- 
body Hamiltonian before diagonalisation, or by transforming the two-body potentials in those partial wave channels 
with deeply-bound forbidden states in a way that preserves phase (spectral) equivalence. We here adopt this second 
approach, and use supersymmetric transformations of the two-body potentials in order to eliminate a required set 
of bound states. All the parameters relative to the method are specified in the namelist b2states and the specific 
characteristics of the two-body bound states to be calculated are defined in b2state (see manual for details). 

Sometimes, it is useful to calculate the two body state generated by a given effective interaction, or explore how to 
adjust the two body interaction to obtain a given binding energy. FaCE allows you to calculate bound states without 
feeding them into the SUSY transformation subroutine (see b2states in the manual for details). 

A. Elimination of Two-body bound states 

A supersymmetric transformation of the set of potentials |18J | enables the removal of an arbitrary bound state, while 
keeping the spectral (^-matrix) equivalence of the initial and the transformed Hamiltonians. In partition k, and in 
each two-body spin-parity channel, the initial Hamiltonian Hq for the interaction of bodies i,j couples N two-body 
channels for total angular momentum jk- The two-body equation is then 



h' 



2 r ^2 e te ,T\^ . J^ 



2fii 



f_ , C(C + 1) 
dr2 



e) K (r) + ^ y„„' (r) 0„, (r) = , (29) 



where n is a channel index set \lx^,^ {si^ Sj)Sxf.] jk}, f = ^ij, ^n = IxkM, 3.nd /^^ is the reduced mass for bodies i,j. 
The threshold energies e„ = eg. +es- are included in the diagonal matrix elements y„„(r), in addition to the couplings 
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defined by Eq. H15() . We start with the real symmetric potential matrix Vo(r) = {Ki'„(r)} at each radius, and repeat 
the following supersymmetric transformation for each bound state p = 1 up to the number of forbidden states P. 

Let the column vector ^p^i{r) = {<l>n{''')} t>e the normalised ground state eigensolution of i/p_i at real energy E\ 
below all thresholds. By applying a double supersymmetric transformation to -ffp-i we obtain a new Hamiltonian Hp 
where the potential matrix V^-i(r) is replaced by Vp{r) 



h' d $^i(r)$^li(r) 



Vp{r) = Vp.,ir) - -- ..4^^\:r .,;„ (30) 



In the case of vanishing coupling between the channels near the origin, it is possible to deduce the behaviour of 
the diagonal parts of Vp{r) at small r. If the diagonal matrix of angular moments {in} has only one lowest element, 
say £1, such that ii < £„ for n = 2,3, ...N, in this channel (index 1) the additional term in the supersymmetric 
transformed potential Vp{r) will have a singularity h {2ii + 3)/(2/iyr^) at small r, which added to the centrifugal 
term ?i ^i(^i + l)/(2/x.yr^) gives a new centrifugal term ?i (£i+2)(^i+3)/(2/x.yr^). The supersymmetric transformation 
Eq. H30|) in this case adds a repulsive core at the origin, by increasing the orbital moment ii by 2 units. In all other 
channels the orbital moments £n are not changed. Physically, this corresponds to the conservation of the oscillator 
quanta A = 2nr + ^ in the system: when reducing the radial quantum number n^ by one unit (removing one level) 
we increase the orbital part £ by two (for the one channel case we satisfy the Levinson's theorem). If the diagonal 
matrix of angular momenta {£n} has several lowest equal elements, the increase of singularity is shared between these 
channels, including their coupling potentials. 

In FaCE, if supersymmetric transformations are used for any partition fc, then the transformations up to Vp{r) 
must be recalculated for all desired two-body channels jk (all jv in namelist b2state). 

V. SOLVING THE FADDEEV EQUATIONS 

The Faddeev equations Eq. (^ are solved, after expanding on the hyper-angular Eq. ^ and hyper-radial Eq. @ 
basis functions, to find square-integrable solutions Eq. (|12|l for eigenenergies E and eigenvectors a. For E < these 
are bound states, whereas for E > the eigen-solutions are 'quasi-bound' states that form a discrete representation 
of the continuum. These quasi-bound wave functions may be used, for example in [10'|, in the calculations of breakup 
as inelastic excitations. 

The physical normalisation of the wave functions Eq. I|12|) is (^|^) = 1. If A^ is the whole normalisation matrix 
A™.Q,., the eigenvectors a are physically normalised when a^Na = 1. Note that some eigenvectors will be found that 

are non-physical, having a^ Na = 0; in these there is a cancellation between different Faddeev components, and they 
must be omitted in all bound or breakup state analyses. 

A. Hamiltonian Reduction procedures 

The complete set Eq. Q of Faddeev equations may be reduced in a number of circumstances. FaCE has the 
option of 'isospin' and 'orthonormal' reductions, which exactly reproduce a physically chosen subset of the eigen- 
solutions, and also 'Feshbach' reduction, which is a method for approximating the effects of high K partial waves on 
the solutions. The choice of the reduction method is made through eqn in the input namelist solve (see manual). 

1. Isospin Reduction 

Suppose bodies j and k are fermions which are isospin states T^ of some particle of isospin Tjk , with Sj = Sk half- 
integral. The requirement of antisymmetrisation under exchange of these bodies is easily satisfied if the partial wave 
set ai only includes those quantum number sets for which l^i + Sxi + Tjk is odd. In this way, wave function components 
that are symmetric under interchange of j and k are eliminated from the basis set for this Faddeev component. 

Furthermore, the remaining Faddeev components '^j and ^fe are isospin mirrors of each other. Just one of these 
wave functions needs to be included explicitly in the equation set to be found numerically, since 

^j = -(-1)^^-^,^*^ (31) 



where Pjk is the operator permuting the coordinates of particles j and k. 
The coupled equations 




V, ^- \ ( ** 

T,+h + V,-E V, vp^. I =0 (32) 



Vk n + h + Vk-E \^ 



k 



are now reduced to 



T, + h + V,-E V^ + V,P,k \ ( *. 

Vu Tk + h + Vk + VkPjk -E )\^k 

The permutation matrix elements are Pjk = ( — 1)('j;j+'S'^j^''2j-s33)^ 



. (33) 



2. Orthonormal Reduction 

Since basis states \i : cti) in Eq. H12|l form an over-complete set, the same set of physical eigen-solutions may be 
found by transforming the basis set into an orthonormal one. 

A rotation matrix C may be found, for example by Gramm-Schmidt orthonormalisation, such that C^ NC ~ I, so 
that the columns of C are vectors that are physically orthonormalised by the norm matrix N. This same rotation 
may be used to transform the Hamiltonian matrix of Eq. I|ll() . Defining D = C^^ = C'^ N, then 

DHCa = H'a' = Ea' 

is an orthogonal transformation of the original eigenvalue problem, with the same eigenvalues. The original eigenso- 
lutions may be regained as a' = Ca. The new matrix H' is real and symmetric; this proves that the eigenvalues of 
Eq. (|11() are real even though H is not symmetric. 

3. Feshbach Reduction 

Another reduction method, also called the semi-adiabatic reduction method, constructs an effective coupling matrix 
at each p value using Feshbach's expression (13 for effective interactions in a subspace. 

Consider the set of N coupled equations for the wave functions X^a k a-s in Eq. IjHJi. Take the subset of the equations 
of this system with largest Ki and core excitation energy e^.^. In these channels, an adiabatic condition might be 
fulfilled, where the hyper-radial kinetic energy Tp is small and can be neglected. Thus we have an option of keeping 
this kinetic energy term in only the subset of channels i = 1, • • • , M, and of neglecting Tp for i~M-\-l,---,N. For 
each p value, let us rewrite the system Eq. ^ in the following matrix form: 

A - £; B \ f x^*"' 

X 



C D - S M y(b) ' ~ ^' ^^^^ 



where A contains the exact Tp + LKi{p) terms, but D contains only the LKi{p) terms. The B and C are the block 
off-diagonal matrices. The solution vectors are x^^'' — (Xi ■ ' ' Xm) and x'"^^ = (xm+i • • ■ Xn)- 
Solving the matrix Eq. H34|l formally we obtain: 

x('') = (D-S)-iCx("), (35) 

and substituting Eq. (|35|l into our system Eq. H34(l we get a reduced subset of coupled equations for x*-"-* 

(A-£;-hB(D-£;)-iC)x^"^ =0. (36) 

^From Eq. (|36l) , we see that the reduction of the coupled equations from A^ x A^ to a smaller M x M set consists 
in adding a 'Feshbach' term B(D — _E)~^C to the effective interaction in the retained subspace. 

Strictly speaking, the Feshbach term should be recalculated for every eigen-energy E, but in practice we calculate 
the Feshbach term once for the fixed 'Feshbach energy' E = Ep, which should be chosen near the eigen-energy of the 
state of most interest, such as the ground state energy. Variables efesh and kmaxf in the input namelist solve are the 
Feshbach energy and the K- value above which the Feshbach approximation is introduced (see manual). 
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B. Diagonalisation procedure 

The subroutine fadco in FaCE evaluates all the potential matrix elements as functions of p. After the above 
possible reductions, the Hamiltonian matrix appearing in the Faddeev Equations Eq. Hll() is determined by hyper- 
radial integrals using the radial basis function Eq. Q. 

For the general eigenvalue solution of Eq. Hll|l . where H is a real matrix not necessarily symmetric, we use the 
subroutine F02AGF from the Nag library, which proceeds via reduction to Hessenberg form, to find all eigensolutions. 

If only selected eigenvalues are required, and the input parameter meigs (introduced in solve) is non-zero and 
less than the dimension of the H matrix of Eq. I)ll|l . a more efficient method is that of inverse iteration. Starting 
from some energy Eq, and some initial guess a''^' for an eigenvector, the solution of the simultaneous linear equations 
{H — Eo)aS'^'> = a*^"~^' for n > 1 will converge to the eigensolution with energy nearest Eq- This method is effective 
for Eq less than the ground state energy, when there are no nearly-degenerate eigenvalues. It may be generalised to 
finding the several eigenvalues nearest Eq by orthogonalising a^") at each iteration to the set of eigenvectors already 
found. 



VI. COMPUTER PROGRAMME AND INPUT MANUAL 

Given the description covered in the previous sections, the FaCE manual is presented as a sequence of namelists 
with explanatory names for the variables. Nevertheless it is useful to remind the user that the parameters delimiting 
the three body space are: the number of Laguerre polynomials Ni, for the hyper-radial part, together with the Jacobi 
polynomials for the hyper-angular part Njac] the maximum angular momentum that are to be taken into account in 
each Faddeev partition lxmax{i)^ lymax{i)^ and the number of X-harmonics KmaxO-, i)- 

The source code is distributed with separate makefiles for Sun f90 compilers (standing alone, or with system Sun 
Performance Library and Nag libraries) and for Linux, where there are Intel ifort and Portland Group pgf90 makefiles, 
the latter optionally with system Lapack and Bias libraries. The suitable one of these makefiles should be renamed 
to 'makefile'. 

FaCE uses F02AGF, MOIDAF and MOIZAF routines from the Nag library, and DGETRF, DGETRS and DGEMM 
from the Bias library. The original source codes for the Nag and Bias subroutines are contained in the package for 
compilation where not otherwise available. 

A. Input namelists 

• &;fname 

nfile[A*20], desc[A*80] 

nfile is the root name of the output files and desc is a heading that should describe and identify the run. 

• &:scale 

amn, he [r*8] 

amn is the mass of the nucleon and he is the planck constant multiplied by the velocity of light in [MeV fm] . 

• &;nuclei 

name(l:3)[A*8], mass(l:3), z(l:3), radius(l:3)[r*8] 

Each of the three interacting nuclei are characterised by their name, mass, charge and radius. 

• &;identical 

id(l :3) [logical], iso(l :3)[r*8] 

If id(j) is true then the interacting pair in partition j are 2 identical particles. This variable affects the choice of 
the basis: if particles are identical, an isospin reduction is performed, iso contains the isospin of the interacting 
pair to be used when id{j) is true, to omit non-antisymmetric partial waves in that partition. 

• &:total 

ngtfint], gtot(l:ngt)[r*8], gparity(l:ngt)[int] 

ngt — + the number of states to be calculated, gtot and gparity hold their total spin and parity(even= +1, 

odd= —1) respectively. 
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• &:particles 

ns(l:3)[int], spin(l:ns,l:3)[r*8], parity ( l:ns,l: 3 )[int], energy(l:ns,l:3)[r*8] 

This namelist defines intrinsic properties of the three bodies. For each of the 3 nuclei: nsQ) specifies the number 
of states to be included for nucleus j. For each of its ns(j) states, one should specify the spin, parity(+l/-l) and 
its energy relative to the energy of the ground state of that nucleus j. 

corek(l:3), def(2:mmuHipoles,l:3)[r*8] 

This namelist contains the electromagnetic information on each of the three bodies. 
corek(j) -^ projection of the spin of nucleus j in its rest frame 
def(q,j) -^ deformation length in the q multipole of nucleus j. 

• &waves 

This namelist contains the definition of the channels per partition. 

sym(3)[A*l] the name (e.g 'T', 'X' or 'Y') of each partition. 

auto (3) [logical] —> if 'T' then FACE generates automatically the channels allowed given maximum quantum 
numbers for each Faddeev partition. Otherwise, quantum numbers for nc{j) channels are explicitly read in. 

If auto is false: 

nc(3)[int] — > number of channels per partition (less than mfchan). 

lx(mfchan,3)[int], ly (mfchan, 3 )[int], lt(mfchan,3) [int], sx(mfchan,3)[r*8], jp(mfchan,3) [r*8] -^ determine 

the quantum numbers associated with all channels for each partition in the following coupling order 

\[{lxJv)lt, {si,Sj)sx\Jp,icy) 

icy(rnjchan,3)[int],icxl(rnfchan,3)[int],icx2(rnfchan,3)[int] -^ the state of the spectator, the first and the second 

interacting particle for the given channel in each partition. The spin of a given state of each of the bodies was 

defined in particles. 

np (mfchan, 3) [int] -^ the number of iiT- harmonics. 

If auto is true: 

lxmax(l:3)[int], lymax(l:3)[int], ltmax(l:3)[int], sxmax(l:3)[r*8], jabmax(l:3)[r*8], 

kmaxa(l :3) [int] , kmax(0:maxl,l:3) [int] — > these establish the maximum quantum numbers allowed in each par- 
tition, kmaxa is the K limit for all partial waves, and may be overridden by particular kmax specified. 

• &poten 

This namelist defines the radial behaviour of the potentials for each interacting pair. 

detail(3)[A *80] — > information on the interaction between the interacting pair in that partition 

typc(3)[A*3], pa(6,3),ps(6,3),pp(6,3),pd(6,3),pf(6,3)fr*8] -^ central interaction 

typso(3)[A*3],pso(6,3),psop(6,3),psod(6,3),psof(6,3)[r*8] -^ spin-orbit interaction 

typss(3) [A*3] ,pss(6,3) ,psss(6,3) ,pssp(6,3) ,pssd(6,3) ,pssf(6,3) -^ spin-spin interaction 

typt(6,3)[A*3], pt(6,3)[r*8] -^ tensor interaction 

rcoul(3) [r*8] ,acoul(3) [r*8] -^ the Coulomb interaction 

lpot(3)[int] — > is useful for /—dependent interactions where there is an ambiguity on the radial form factor that 

should be used for off diagonal couplings. If Ipot — 0, the radial form factor corresponding to the minimum li, If 

is used; Ipot = 1, the average is taken; Ipot — 2 the maximum is used; and lpot=3 the final If (corresponding to 

the left hand side of the matrix element) is used. Finally, when Ipot > 10, the radial form factor for off-diagonal 

coupling is determined by I — Ipot — 10, throughout the whole calculation, leaving the monopolc terms untouched. 

The form factor for the potentials between the interacting pair in each partition is specified by type (gau, ws, 
rnp, rnn, nul) and the potential parameters for each partial wave (s,p,d,f and a or no extra letter for all). If 
the type is gau then the interaction is the sum of 3 Gaussians: 



Vgauir) == Yl Pa(fc,i)exp 



pa{k + 1, i) 



fe=l,3,5 

If type is ws then the interaction is the sum of 2 Woods-Saxon: 

' r — pa{k + 1, i) 



VLir)^ Ep«(^'*) 



fc=l,4 



1 + exp 



pa{k + 2, i) 



(37) 



(38) 
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For the spin-orbit interaction, if type is 'ws' then the form factor is given by the derivative of two Woods-Saxon: 



VL{r) 



pso(k,i) 



exp( 



r— pso(fc-f-l,i) 
pso{k'\-2,i) 



^ r psoffc -I- 2, i) fi +pxr)('^-^£^^^^-ti4i 

k=lA ^ ^ ' ' V^ ^ '^^Vy pso(k+2,i} 



(39) 



The Coulomb interaction for the interacting pair in partition (i) is that of a uniform sphere with radius rcoul(i) 
and diffuseness acoul(i), screened at radius rscreen with a Fermi function of diffuseness ascreen. 

The operator for the tensor force is 12512 as defined by Brink and Satchler ^^. The operator for the spin-spin 
force is the dot product of the spins of the interacting pair Sj ■ Sk . The operator for the spin-orbit is defined in 
gamso. If the deformation of one of the interacting particles in non zero then higher order multipoles will be 
automatically added to the monopole interaction based on a spherical harmonic decomposition of a deformed 
field. 

• &:pot3b 

typ3b,s3b(ngt) ,r3h(ngt) ,a3h(ngt) ,gtvary — > specifies the parameters for the diagonal 3-body potential V3(p) if 
typ3b 7^ nul. If gtvary, then ijt = 1 below, otherwise ijt is the J^ index 1 ... ngt. If typSb ~ gau, then 



Vsip) = s36(ijt)exp 



r3b{ijt) 



(40) 



If typSb — ws, then 



Vsip) = s3b{ijt) 



1 + exp 



p - r3b{ijt) \ 
a3b{ijt) J 



(41) 



• &:gamso 

gamso 1, gamso2 —>■ the spin-orbit matrix elements are calculated using the following operator Til^ ■ si - 
for each partition 



'T2lx-S2 



• &grids 

rr [r*8], nlag,njac [int] 

This section contains radii for the expansions used. 

rr — > scaling parameter po for the Laguerre basis 

nlag -^ number of Laguerre quadrature points for the p coordinate. 

njac -^ number of Jacobi polynomials 



• &;trace 

pripotjVadia [logical] 

Printing options: pripot prints the potential matrix elements, vadia prints the diagonalised coupling eigenvalues 

(energy surfaces). All of these are printed in the output file with extension lis. 



&:b2states 

n2states [int] , dx,xmax, [r*8], ipc,lmax,nk [int], rnode,de [r*8], 

Find n2states two-body states in the two-body potentials. Use radial grid to xmax with steps dx. 

Imax, nk and mode, de are default values for each b2state namelist below. 



The ipc. 



• &:b2state (repeated n2states times) 
pair, kind [int], de [r*8] ipc [int], test [logical], n,nvchan,l,lmax [int], s,jv, mode [r*8], search, rescale [log- 
ical], eigen, potential, fermi [r*8], nomit [int],omit_ll:nomit [int] ,omit_sl:nomit,omit_jl:nomit [r*8], 
omit_cl 1 :nomit, omit_c21 :nomit [int], 

Find two-body eigenstate in the potential pair, of fcmd='occup' to be used for pauli blocking via the susy 
transformation; A;md='transfer' if one needs to check the properties of a particular two body state, or fcmd^'pot' 
if one needs to calculate numerically the potential for a particular two body partial wave set, without excluding 
it. This last option is needed, because whenever there are any occupied states, the potentials for all partial 
waves sets in that partition need to be calculated numerically. 
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ipc = trace level, test=T to ignore this state after finding it. 

Wave functions will have n nodes in channel £—1 up to radius mode, from a set of £ <lmax using coupling order 

Ki {siS2)s; jv). The eigenenergy is eigen in monopole potential multiplied by potential, where search=^W or 'V 

to search for energy or potential factor respectively. Energies eigen are negative for bound states. Only bound 

states can be found when search ='E'. Use nomit>0 to specifically omit some partial waves from the coupled 

channels set. 

fermi <0, to exclude bound states up to that valence energy, 

fermi >0, to exclude the nint (/ermi) number of lowest-energy bound states. 

&solve 

eig, eimin,eimax,efesh [r*8], eqn[A*l], cfiles [logical], nbmax, meigs(l:ngt),kmaxf [int] 

This namelist is related with the type of equation to be solved: 

eig, eimin and eimax are the target, minimum and maximum eigenenergies to search for states, 

eqn asks for the reduction of the full equation: eqn= 'F' stands for Faddeev, eqn=sym(i) (defined in waves) 

perforins an orthonormal transformation to the i basis, 

nbmax = number of functions in the radial expansion, must be < nlag, 

meigs is the number of eigenstates to calculate {meigs=0 is to find all eigenstates). 

kmaxf > for Feshbach reduction of coupled equations at each hyper-radius to K < kmaxf, using eigenenergy 

estimate efesh. 

cfiles = 'T', to write mel and spec files to be fed into an independent program of the coupled equations (e.g 

the program sturmxx 20]). 



B. Outputs 

• standard output 

The standard output contains the information about the three nuclei, the partial waves to be included, the 
two body potentials, the parameters used in the expansions. If the run uses supersymmetric potentials, the 
details regarding the two body bound states to be excluded are printed out. Next, FaCE prints the angular 
momentum information about all the possible channels, the Gauss-Laguerre grid, the details about the reductions 
performed and the corresponding new reduced set of channels. Finally the energy, the radii, and the probabilities 
associated with each channel are given for each calculated state: values for L-summed probabilities, and summed 
probabilities for each core-state are also included. As JJ coupling is easier to compare with the shell model 
basis, FaCE performed the LS-JJ transformation and prints out the probabilities of the main JJ components at 
the very end of the file. 

The following files are produced with filename = nfile in the fname namelist. 

• filename. viff 

This file contains the hyper-radial wavefunction for the states J'^ calculated. It first contains the channels 
that are included in this output (very small components are left out) followed by the wavefunction in format 
r, wf{i, r). This file can be easily plotted. 

• filename. nl 

FaCE rewrites into filename.nl the input as is read. 

• filename. lis 

This file contains extensive information on the various steps of the calculation. It contains the two body 
potentials for each partition, the algebra matrix elements presented in section III, the various radial potential 
couplings for the various hyper-radii belonging to the Gauss-Laguerre grid, the normalization and permutation 
matrices, and the probability of the various configurations in long format. 
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VII. EXAMPLES OF CALCULATIONS 
A. i^Be 

Input file hel2gptdefk4-in and a shortened version of the output file bel2gptdefk4-out are provided below. The full 
files are included in the electronic file distribution. This example models ^^Be as a three-body cluster of two neutrons 
outside a ^°Be core. The core is deformed and allowed to excited to its first 2+ state. This example is similar to that 
in 1^ although here we use shallow core-n potentials, to most simply avoid Pauli-forbidden two-body states. 

1. Example: bel2gptdefk4-in 

kfname nf ile='bel2gptdefk4' desc= 'bel2gptdefk4: n+n+belO using gptnn and belO-n , Kmax=4' / 

fescale amn=939. hc=197.3/ 

fenuclei name= 'n' , 'n' , 'lObe' mass= 1 1 10 z= 4 radius=0 2.30 / 

feldentical id=F, F, F, lso=0.5, 0.5, 0. / 

fetotal ngt = 1, gtot(l)=0.0, gparity(l)=+l / 

&p articles 

ns(l)=l, spin(l,l)= 0.5, parity(l,l)=l, energyCl , 1)=0 .0 , 

ns(2)=l, spin(l,2)= 0.5, parity (1,2)=1 , energyCl ,2)=0 .0 , 

ns(3)=2, spin(l,3)= 0.0, parity (1,3)=1 , energyCl ,3)=0.0, 

spinC2,3)= 2.0, parity C2, 3) =1, energyC2,3)=3.368/ 
&em 

corekCl)=0.5 defC2,l)=0.0 qniomCl)=0.0 MmomCl)=0.0 

corekC2)=0.5 defC2,2)=0.0 ClmomC2)=0.0 MmomC2)=0.0 

corekC3)=0.0 def C2,3)=1.6638 qmom(3)=0.0 Mmom(3)=0.0 defC4,3)=0 / 

fewaves autoCl)=T, kmaxaCl)=4, lxmaxCl)=2, 

autoC2)=T, kniaxaC2)=4, lxmaxC2)=2, 
autoC3)=T, kniaxaC3)=4, lxmaxC3)=2/ 
fepoten 

detailCl) ='n+10be' typcCl) ='ws' 

psCl,l)=-10.14 psC2,l)=2.736, psC3,l)=0.67 
ppCl,l)=-24.24 ppC2,l)=2.736, ppC3,l)=0.67 
pdCl,l)=-10.14 pdC2,l)=2.736, pdC3,l)=0.67 
lpotCl)=0 
typsoCl)='ws' 

psopCl,l)=+25.72 psopC2,l)=2.736, psop(3, 1)=0. 67 
psodCl,l)=-25.72 psodC2,l)=2.736, psodC3, 1)=0. 67 
typssCl)='nul' typtCl) ='nul' 

detailC2) ='10be+n' typcC2) ='ws' 

psCl,2)=-10.14 psC2,2)=2.736, psC3,2)=0.67 

ppCl,2)=-24.24 ppC2,2)=2.736, ppC3,2)=0.67 

pdCl,2)=-10.14 pdC2,2)=2.736, pdC3,2)=0.67 

lpotC2)=0 
typsoC2)='ws' 

psopCl,2)=+25.72 psopC2,2)=2.736, psop(3,2)=0.67 

psodCl,2)=-25.72 psodC2,2)=2.736, psod(3,2)=0.67 
typssC2)='nul' typtC2) ='nul' 

detailC3) = 'gptnn' typcC3) ='gau' 

psCl,3)= 560.0 psC2,3)=0.8109 psC3,3)=-390 .7 

psC4,3)=1.031 psC5,3)=-1.501 psC6,3)=3. 205 

ppCl,3)= 9.335 ppC2,3)= 1.184 ppC3,3)= -1.37 

ppC4,3)= 2.099 ppC5,3)= 0.1663 ppC6,3)= 3.193 

pdCl,3)= 560.0 pdC2,3)= 0.8109 pdC3,3)= -390.7 
pdC4,3)= 1.031 pdC5,3)= -1.501 pdC6,3)= 3.205 
typsoC3)='gau' 

psoCl,3)= -114.5 psoC2,3)=0.9296 
typssC3)='nul' typtC3) ='gau' 

ptCl,3)= 12.24 ptC2,3)= 1.539 ptC3,3)= -31.64 

ptC4,3)= 0.4039 ptC5,3)= 0.8111 ptC6,3)= 3.015 
/ 
&pot3b typ3b='nul', s3b=0., r3b=3.9 / 

gamsol=l,0,l, gamso2=0 , 1 , 1 / 
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fegrids rr=0.3 nlag=20 njac=40 / Methods: 

fetrace pripot='T' / 

&b2states n2states=0 dx=0.002 xmax=15. / 

&solve eimin = -5.0, eimax=3, eqn='F' nbmax=10 meigs=l / 



2. Output: bel2gptdefk4-out 



FACE: version 0.12e 

Case file: bel2gptdefk4 

bel2gptdefk4: n+n+belO using gptnn and belO-n , Kmax=4 



with constants: 

unit mass = 939.00 MeV; he 



197.300 MeV.fm => h2sm = 20.7281 



12 3 
Nuclei : n n lObe 

masses 1.0000 1.0000 10.0000 
charges 0.0 0.0 4.0 



radii 



0.0000 0.0000 2.3000 



Identical 23: F, 31: F, 12: F, 
isospin 0.5 0.5 0.0 

1 coupled channels J, pi sets: 0.0+ # 1 



Particle 1: n has 1 states: 

— spin, parity, energy = 0.5+ 0.0000 MeV 

— Intrinsic K = 0.5 Quadrupole moment = 

Particle 2: n has 1 states: 

— spin, parity, energy = 0.5+ 0.0000 MeV 

— Intrinsic K = 0.5 Quadrupole moment = 



0.000 Magnetic moment = 0.000 



0.000 Magnetic moment = 0.000 



Particle 3: lObe has 2 states: 

— spin, parity, energy = 0.0+ 0.0000 MeV 2.0+0 3.3680 MeV 

— Intrinsic K = 0.0 Quadrupole moment = 0.000 Magnetic moment = 0.000 

— deformation lengths = 1.66380 



Partial waves 
Component 1 X 
Component 2 Y 
Component 3 T 



lx,ly,lt <= 2 10 10, sx,jp <= 2.510.0, Kmax(all,0: Ix) = 4-1-1 -1 
lx,ly,lt <= 2 10 10, sx,jp <= 2.510.0, Kmax(all,0: Ix) = 4-1-1 -1 
lx,ly,lt <= 2 10 10, sx,jp <= 1.010.0, Kmax(all,0: Ix) = 4-1-1 -1 



POTENTIALS 

Potential 1 between n and lObe : 
n+lObe 

Central potential of type 'ws ' , 
for s-waves: -10.14000 2.73600 0.67000 0.00000 0.00000 0.00000 
for p-waves: -24.24000 2.73600 0.67000 0.00000 0.00000 0.00000 
for d-waves: -10.14000 2.73600 0.67000 0.00000 0.00000 0.00000 
using Ipot = 

Spin-orbit potential of type 'ws ' , 
for p-waves: 25.72000 2.73600 0.67000 0.00000 0.00000 0.00000 
for d-waves: -25.72000 2.73600 0.67000 0.00000 0.00000 0.00000 
acting on n with factor 1.0000, on lObe with factor 0.0000 

Spin-spin potential of type 'nul' , 

Tensor potential of type 'nul' , 



Potential 2 between lObe 
lObe+n 
Central potential of type 
for s-waves: -10.14000 



and n 



ws ' , 
2.73600 



0.67000 0.00000 0.00000 0.00000 
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for p-waves: -24.24000 2.73600 0.67000 
for d-waves: -10.14000 2.73600 0.67000 
using Ipot = 
Spin-orbit potential of type 'ws ' , 
for p-waves: 25.72000 2.73600 0.67000 
for d-waves: -25.72000 2.73600 0.67000 
acting on lObe with factor 0.0000, on n 
Spin-spin potential of type 'nul' , 
Tensor potential of type 'nul' , 



0.00000 
. 00000 



0.00000 
0.00000 



0.00000 
0.00000 



0.00000 
0.00000 



with factor 



0.00000 
0.00000 



0.00000 
0.00000 
1.0000 



Potential 3 between n and n 

gptnn 
Central potential of type 

for s-waves: 560.00000 

for p-waves: 9.33500 

for d-waves: 560.00000 

using Ipot = 
Spin-orbit potential of type 'gau' 

for all-waves: -114. 50000 0.92960 
acting on n with factor 1 
Spin-spin potential of type 'nul' 
Tensor potential of type 'gau' 



gau' , 
0.81090-390.70000 
1.18400 -1.37000 
0.81090-390.70000 



0.00000 
0000, on n 



1.03100 
2.09900 
1.03100 



-1.50100 

0.16630 

-1.50100 



0.00000 0.00000 
with factor 



20500 
19300 
20500 



0.00000 
1.0000 



12.24000 



1 . 53900 



-31.64000 



0.40390 



0.81110 



3.01500 



Three-body potential of type 'nul' 



METHOD parameters: 

Hyperradial parameters = 0.3000 20 
Hyperangular points = 40 

nbmax = 10 , eig = -5.000000 

nlag = 20 recommend: >> nbmax = 10 

njac = 40 recommend: >> kmax/2 = 2 

(but much more, for repulsive-core interactions) 

Equations to solve = F=F (F=Faddeev, T,X, Y-Schrodinger+orthoNormal) 

For 0.0+ search for 1 eigenstates near -5.000 MeV, 
examine those between -5.000 k 3.000 MeV 

* Coupled channels set 1 for J, pi = 0.0+ * 



Faddeev channel numbers required: 



30 



30 



30 



X: 

ig 1 1 


2 


3 


4 5 6 


7 8 


9 


10 


11 


12 


13 


14 


15 16 


i 1 1 


2 


3 


4 5 6 


7 8 


9 


10 


11 


12 


13 


14 


15 16 


K 1 


2 


2 


2 4 4 


4 4 


4 


2 


2 


2 


2 


2 


2 2 


LI 





1 


1 


1 





2 


1 


2 


2 


2 


2 2 


sx 1 0.5 


0.5 


0.5 


0.5 0.5 0.5 


0.5 0.5 


0.5 


1.5 


1.5 


1.5 


1.5 


2.5 


2.5 2.5 


Ix 1 


1 


1 


2 2 


1 1 








1 


1 


2 





1 2 


ly 1 


1 


1 


2 2 


1 1 





2 


1 


1 





2 


1 


JP 1 0.5 


0.5 


0.5 


0.5 0.5 0.5 


0.5 0.5 


0.5 


0.5 


0.5 


0.5 


0.5 


0.5 


0.5 0.5 


iz 1 1 


1 


1 


1 1 1 


1 1 


1 


2 


2 


2 


2 


2 


2 2 



« 



54 lines deleted 



» 



All Gauss-Laguerre points: 

radii: 0.252 0.500 0.813 

3.452 4.217 5.073 

11.036 12.680 14.549 

Gauss-Laguerre points for nbmax = 

radii: 0.451 0.901 1.478 

6.952 8.896 11.527 



1.194 


1.646 


2.170 


2.771 


6.026 


7.086 


8.263 


9.573 


16.712 


19.307 


22.679 




10 : 








2.196 


3.071 


4.127 


5.401 
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Calculating HH up to n = 
Calculating matrix elements 
do potentials for partition 
do potentials for partition 
do potentials for partition 
allocate ww section of 
allocate ww files of 2 * 



11,12 = 



1 
2 
3 

1 Mb 

2 Mb 

1 . 048809 

1 . 048809 

1.414214 



For basis 1 , dnn 

For basis 2 , dnn 

For basis 3 , dnn 

eqn,id(:)=F F F F 

allocate wm array of 8 Mb 

Matrix elements found in Gauss-Laguerre Basis 
Diagonalise matrix of size 990 



Search for 




1 


eigenstates near -5.000000 


Inverse iteration to 


find 


1 eigenstates nearest -5.00 


iteration 


1 gives 


-4 


70616, change = -4.71E+00 


iteration 


2 gives 


-4 


87763, change = -1.71E-01 


iteration 


3 gives 


-4 


87759, change = 3.51E-05 


iteration 


4 gives 


-4 


87759, change = 6.10E-07 



Found 0.0+ Evals at: 



-4.87759 



Bound 0.0+ eig 1 = 
rms <rho> = 4.514 fm. 



-4.877592 

rms matter radius = 



2.471 fm. 



Probability norms in eigenstate: 



K L sx Ix ly jp i# 



X: 



0.5 

0.5 

1 0.5 
0.5 
0.5 



0. 

1 0. 

1 0. 
0. 

2 0. 



25 lines deleted 



1 
1 
1 
1 
1 
» 



normalised 

0.06145327 
0.04085356 
0.01151844 
0.00336356 
. 00004349 



permuted 

0.58052013 
0.29732071 
. 04324449 
0.00247767 
0.00563919 



P(S) = 0.583931, P(S') = 0.303015, P(P) = 0.045894, P(D) = 0.067025 



Y: 
31 
32 
33 
34 
35 



0.5 

0.5 

1 0.5 
0.5 
0.5 



0. 

1 0. 

1 0. 
0. 

2 0. 



25 lines deleted 



1 
1 
1 
1 
1 
» 



0.06145327 
0.04085356 
0.01151844 
0.00336356 
. 00004349 



0.58052013 
0.29732071 
. 04324449 
0.00247767 
0.00563919 



P(S) = 0.583931, P(S') = 0.303015, P(P) = 0.045894, P(D) = 0.067025 



61 





0.0 





0.0 1 




0.07138963 


0.58052013 


62 


2 


0.0 


1 


1 0.0 1 




0.00000000 


0.00000000 


63 


2 


0.0 





0.0 1 




0.02395303 


0.29979838 


64 


2 


1 1.0 


1 


1 0.0 1 




0.00005476 


. 04324449 


65 


4 
. 2 


0.0 
5 lines 


2 
de 


2 0.0 1 
leted . . 


. » 


0.00000006 


0.00202213 



P(S) = 0.884924, P(S') = 0.002022, P(P) = 0.045894, P(D) = 0.067137 

norm of X channels = 0.999866 
norm of Y channels = 0.999866 
norm of T channels = 0.999977 



J J coupling: X 

[01/2+0 0/2 ]0.5, 0.5 

[11/2+1 2/2 ]0.5, 0.5 

[13/2+1 2/2 ]0.5, 0.5 

[11/2+1 2/2 ]0.5, 0.5 



0.0 
0.0 
0.0 
0.0 



0.58393105 
0.23452232 
0.10610143 
0.02629508 



18 



[13/2+1 2/2 ]0.5, 0.5; 0.0 : 0.03494796 
ProbClObe in state 1) = 0.930569 from Jab 0.000000 
ProbClObe in state 2) = 0.069297 from Jab 0.000000 
Total norm in jj basis = 0.999866 



0.930569 0.000000 
0.069297 0.000000 



J J coupling: Y 










[00/2+0 1/2 ]0.5, 





5; 


0.0 


0.58393105 


[12/2+1 1/2 ]0.5, 





5; 


0.0 


0.02139192 


[12/2+1 3/2 ]0.5, 





5; 


0.0 


0.31923183 


[12/2+1 1/2 ]0.5, 





5; 


0.0 


0.01520326 


[12/2+1 3/2 ]0.5, 





5; 


0.0 


0.04597204 


ProbClObe in state 1) 


= 


0. 


930569 


from Jab 0.000000 


ProbClObe in state 2) 


= 


0. 


069297 


from Jab 0.000000 


Total norm in jj basis 


= 


C 


.999866 


J J coupling: T 










[01/2+0 1/2 ]0.0, 





0; 


0.0 


0.88492424 


[11/2+1 1/2 ]0.0, 





0; 


0.0 


0.02908190 


[13/2+1 3/2 ]0.0, 





0; 


0.0 


0.01454095 


[01/2+2 3/2 ]2.0, 


2 


0; 


0.0 


0.02071231 


[01/2+2 5/2 ]2.0, 


2 


0; 


0.0 


0.03106845 


ProbClObe in state 1) 


= 


0. 


930569 


from Jab 0.930569 


ProbClObe in state 2) 


= 


0. 


069408 


from Jab 0.000000 


Total norm in jj basis 


= 


C 


.99997" 




round state 0.0+ ENGY 


= 




-4.877E 


)92 



0.930569 0.000000 
0.069297 0.000000 



0.000000 0.000000 
0.000000 0.069408 



B. "He 

Input file he6psk06.in and complete output file he6psk06.out are provided in the electronic file distribution. This 
models two neutrons outside an inert '^He core with a Gaussian-shape n-^He potential, with SUSY elimination of the 
Os bound state in the s-wave potential, and a GPT nn potential, as in pi3lll7|. 



C. «B 



Input file h8griglk6.in and complete output file h8griglk6.out are provided in the electronic file distribution. This 
models ®B as a three-body cluster of -^116, "^He and a proton, as in 
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